    module recon_mod
      use spark
      use constants_mod
      use params_mod
      use weno_coef_mod
      use math_mod
      implicit none
      
      real   (r_kind), dimension(:), allocatable :: quad_wts_1d
      real   (r_kind), dimension(:), allocatable :: quad_wts_2d
      
      integer(i_kind) :: maxRecCells
      integer(i_kind) :: maxRecTerms
  
      real   (r_kind) :: recCoef
      real   (r_kind) :: recdx
      real   (r_kind) :: recdy
      real   (r_kind) :: recdV
      
      ! For WENO
      ! WENO stencil type, record the indices of cells in corner
      integer(i_kind) :: nWenoStencil
      integer(i_kind) :: nWenoCells   ! nCells in each WENO reconstruction stencil
      integer(i_kind) :: nWenoTerms   ! nTerms in each WENO reconstruction stencil
      integer(i_kind) :: nWenoType   
      integer(i_kind) :: nWenoPoints  ! nPoints on cell those are reconstructed by WENO
      integer(i_kind) :: nWenoLack    ! number of lack cells
      
      integer(i_kind), dimension(:,:,:), allocatable :: wenoType
      
      real   (r_kind), dimension(:,:,:), allocatable :: wenoCoef
      logical        , dimension(:,:  ), allocatable :: availableStencil
      
      real   (r_kind), dimension(:,:,:), allocatable :: AWENO    ! (nWenoStencil,nWenoTerms,nWenoTerms)
      real   (r_kind), dimension(:,:,:), allocatable :: iAWENO   ! (nWenoStencil,nWenoTerms,nWenoTerms)
      integer(i_kind), dimension(:,:  ), allocatable :: wenoIdx  ! (nWenoStencil,nWenoCells)
      
      real   (r_kind), dimension(:,:,:), allocatable :: iAPoly   ! (nWenoType,nWenoPoints,nWenoTerms)
      real   (r_kind), dimension(:,:,:), allocatable :: iAPoly_x ! (nWenoType,nWenoPoints,nWenoTerms)
      real   (r_kind), dimension(:,:,:), allocatable :: iAPoly_y ! (nWenoType,nWenoPoints,nWenoTerms)
      integer(i_kind), dimension(:,:  ), allocatable :: existPolyTerm
      integer(i_kind), dimension(:,:  ), allocatable :: wenoLack
      
      integer(i_kind), dimension(:,:), allocatable :: xDirWENO
      integer(i_kind), dimension(:,:), allocatable :: yDirWENO
      
      integer(i_kind), dimension(:,:,:), allocatable :: wenoLIdx
      integer(i_kind), dimension(:,:,:), allocatable :: wenoRIdx
      integer(i_kind), dimension(:,:,:), allocatable :: wenoBIdx
      integer(i_kind), dimension(:,:,:), allocatable :: wenoTIdx
      integer(i_kind), dimension(:,:,:), allocatable :: wenoQIdx
      
      integer(i_kind) :: stencil_width_sub
      
      real(r_kind), dimension(:  ), allocatable :: xw
      real(r_kind), dimension(:  ), allocatable :: yw
      real(r_kind), dimension(:,:), allocatable :: ps ! position polynomial of reconstructed points
      real(r_kind), dimension(:,:), allocatable :: px ! position polynomial of reconstructed points for x derivative 
      real(r_kind), dimension(:,:), allocatable :: py ! position polynomial of reconstructed points for y derivative 
      ! For WENO
      
    contains
      subroutine recon_init(mesh)
        type(cubed_sphere_mesh_type), intent(in), target :: mesh
        
        integer(i_kind) :: i,j,k,iStencil,iCOS
        integer(i_kind) :: iQP,jQP,iPOC,iPOE,jPOE
        
        integer(i_kind) :: xdir,ydir
        
        real(r_kind), dimension(:,:), allocatable :: quad_wts_tmp_1d
        real(r_kind), dimension(:,:), allocatable :: quad_wts_tmp_2d
        
        allocate(quad_wts_1d(nPointsOnEdge   ))
        allocate(quad_wts_2d(nPointsOnEdge**2))
        
        allocate(quad_wts_tmp_1d(nPointsOnEdge,1            ))
        allocate(quad_wts_tmp_2d(nPointsOnEdge,nPointsOnEdge))
        
        call mesh%get_params(weq=quad_wts_1d)
        
        quad_wts_tmp_1d(:,1) = quad_wts_1d
        
        quad_wts_tmp_2d = matmul(quad_wts_tmp_1d,transpose(quad_wts_tmp_1d))
        
        k = 0
        do j = 1,nPointsOnEdge
          do i = 1,nPointsOnEdge
            k = k + 1
            quad_wts_2d(k) = quad_wts_tmp_2d(i,j)
          enddo
        enddo
        
        maxRecCells = stencil_width**2
        maxRecTerms = maxRecCells
    
        recCoef = 1.
        recdx   = 1. / ( dx * recCoef )
        recdy   = 1. / ( dy * recCoef )
        recdV   = 1. / ( recCoef**2 )
        
        if(trim(recon_h_scheme)=='WENO2D')then
          if(stencil_width==3)then
            nWenoStencil = 4
            nWenoCells   = 4
          else
            nWenoStencil = ( stencil_width - 2 )**2
            nWenoCells   = 9
          endif
          
          nWenoTerms  = nWenoCells
          nWenoType   = recBdy**2 + 1
          nWenoPoints = 4 * nPointsOnEdge + nPointsOnEdge**2
          nWenoLack   = recBdy**2
          
          allocate( wenoCoef        (nWenoType,nWenoPoints,nWenoStencil) )
          allocate( availableStencil(nWenoType,nWenoStencil) )
          allocate( AWENO           (nWenoStencil,nWenoTerms ,nWenoTerms) )
          allocate( iAWENO          (nWenoStencil,nWenoTerms ,nWenoTerms) )
          allocate( iAPoly          (nWenoType,nWenoPoints,maxRecCells) )
          allocate( iAPoly_x        (nWenoType,nWenoPoints,maxRecCells) )
          allocate( iAPoly_y        (nWenoType,nWenoPoints,maxRecCells) )
          allocate( existPolyTerm   (nWenoType,maxRecCells) )
          allocate( wenoIdx         (nWenoStencil,nWenoCells) )
          allocate( wenoLack        (nWenoType,nWenoLack) )
          allocate( xw              (nWenoPoints) )
          allocate( yw              (nWenoPoints) )
          allocate( ps              (nWenoPoints,nWenoTerms) )
          allocate( px              (nWenoPoints,nWenoTerms) )
          allocate( py              (nWenoPoints,nWenoTerms) )
          
          allocate( wenoLIdx(nPointsOnEdge    , -1:1, -1:1) )
          allocate( wenoRIdx(nPointsOnEdge    , -1:1, -1:1) )
          allocate( wenoBIdx(nPointsOnEdge    , -1:1, -1:1) )
          allocate( wenoTIdx(nPointsOnEdge    , -1:1, -1:1) )
          allocate( wenoQIdx(nQuadPointsOnCell, -1:1, -1:1) )
          
          call calc_weno_coef(wenoCoef,availableStencil,iAWENO,iAPoly,iAPoly_x,iAPoly_y,existPolyTerm,wenoIdx,xw,yw,stencil_width_sub,wenoLack,stencil_width,nPointsOnEdge)
          
          call calc_rectangle_poly_matrix      (stencil_width_sub,stencil_width_sub,nWenoPoints,xw,yw,ps)
          call calc_rectangle_poly_deriv_matrix(stencil_width_sub,stencil_width_sub,nWenoPoints,xw,yw,px,py)
          
          ! Calculate rematch positon for quadrature points in each corner
          xdir = 1
          ydir = 1
          do iPOE = 1,nPointsOnEdge
            wenoLIdx(iPOE,xdir,ydir) =  nPointsOnEdge*0 + iPOE
            wenoRIdx(iPOE,xdir,ydir) =  nPointsOnEdge*1 + iPOE
            wenoBIdx(iPOE,xdir,ydir) =  nPointsOnEdge*2 + iPOE
            wenoTIdx(iPOE,xdir,ydir) =  nPointsOnEdge*3 + iPOE
          enddo
          
          iPOE = 0
          do jQP = 1,nPointsOnEdge
            do iQP = 1,nPointsOnEdge
              iPOE = iPOE + 1
              wenoQIdx(iPOE,xdir,ydir) =  nPointsOnEdge*4 + iPOE
            enddo
          enddo
          
          xdir =-1
          ydir = 1
          do iPOE = 1,nPointsOnEdge
            wenoLIdx(iPOE,xdir,ydir) =  nPointsOnEdge*1 + iPOE
            wenoRIdx(iPOE,xdir,ydir) =  nPointsOnEdge*0 + iPOE
            wenoBIdx(iPOE,xdir,ydir) =  nPointsOnEdge*3 - iPOE + 1
            wenoTIdx(iPOE,xdir,ydir) =  nPointsOnEdge*4 - iPOE + 1
          enddo
          
          iPOE = 0
          do jQP = 1,nPointsOnEdge
            do iQP = 1,nPointsOnEdge
              iPOE = iPOE + 1
              wenoQIdx(iPOE,xdir,ydir) =  nPointsOnEdge*4 + jQP * nPointsOnEdge + xdir * ( iQP -1 )
            enddo
          enddo
          
          xdir = 1
          ydir =-1
          do iPOE = 1,nPointsOnEdge
            wenoLIdx(iPOE,xdir,ydir) =  nPointsOnEdge*1 - iPOE + 1
            wenoRIdx(iPOE,xdir,ydir) =  nPointsOnEdge*2 - iPOE + 1
            wenoBIdx(iPOE,xdir,ydir) =  nPointsOnEdge*3 + iPOE
            wenoTIdx(iPOE,xdir,ydir) =  nPointsOnEdge*2 + iPOE
          enddo
          
          iPOE = 0
          do jQP = 1,nPointsOnEdge
            do iQP = 1,nPointsOnEdge
              iPOE = iPOE + 1
              wenoQIdx(iPOE,xdir,ydir) =  nWenoPoints + ydir * jQP * nPointsOnEdge + xdir * iQP
            enddo
          enddo
          
          xdir =-1
          ydir =-1
          do iPOE = 1,nPointsOnEdge
            wenoLIdx(iPOE,xdir,ydir) =  nPointsOnEdge*2 - iPOE + 1
            wenoRIdx(iPOE,xdir,ydir) =  nPointsOnEdge*1 - iPOE + 1
            wenoBIdx(iPOE,xdir,ydir) =  nPointsOnEdge*4 - iPOE + 1
            wenoTIdx(iPOE,xdir,ydir) =  nPointsOnEdge*3 - iPOE + 1
          enddo
          
          iPOE = 0
          do jQP = 1,nPointsOnEdge
            do iQP = 1,nPointsOnEdge
              iPOE = iPOE + 1
              wenoQIdx(iPOE,xdir,ydir) =  nWenoPoints + ydir * ( jQP - 1 ) * nPointsOnEdge + xdir * ( iQP - 1 )
            enddo
          enddo
        endif
        
      end subroutine recon_init
      
      function edge_quadrature(q)
        real(r_kind) :: edge_quadrature
        real(r_kind) :: q(nPointsOnEdge)
        
        edge_quadrature = dot_product(quad_wts_1d,q)
        
      end function edge_quadrature
      
      function cell_quadrature(q)
        real(r_kind) :: cell_quadrature
        real(r_kind) :: q(nPointsOnEdge**2)
        
        cell_quadrature = dot_product(quad_wts_2d,q)
        
      end function cell_quadrature
    
      subroutine WENO(q,wenoType,qrec,dqdx,dqdy)
        real   (r_kind), intent(in )           :: q   (maxRecCells)
        integer(i_kind), intent(in )           :: wenoType
        real   (r_kind), intent(out), optional :: qrec(nWenoPoints)
        real   (r_kind), intent(out), optional :: dqdx(nWenoPoints)
        real   (r_kind), intent(out), optional :: dqdy(nWenoPoints)
        
        real   (r_kind), parameter :: eps   = 1.e-15
        real   (r_kind), parameter :: theta = 3
        
        integer(i_kind) :: iStencil,jStencil,iCell,iPoint,iTerm,iCount
        
        real(r_kind), dimension(nWenoCells ) :: qC
        
        real(r_kind), dimension(nWenoStencil,nWenoTerms ) :: a  ! coef of reconstruction polynomial
        real(r_kind), dimension(nWenoStencil            ) :: SI
        real(r_kind), dimension(nWenoStencil,nWenoPoints) :: w
        real(r_kind), dimension(nWenoPoints ,nWenoTerms ) :: wa ! weighted polynomial coef
        
        real(r_kind), dimension(nWenoStencil) :: rp
        real(r_kind), dimension(nWenoStencil) :: rn
        
        real(r_kind), dimension(nWenoStencil,nWenoPoints) :: wp
        real(r_kind), dimension(nWenoStencil,nWenoPoints) :: wn
        
        real(r_kind), dimension(nWenoPoints ,nWenoTerms ) :: wap ! weighted polynomial coef
        real(r_kind), dimension(nWenoPoints ,nWenoTerms ) :: wan ! weighted polynomial coef
        
        real(r_kind), dimension(nWenoPoints) :: qrecp
        real(r_kind), dimension(nWenoPoints) :: qrecn
        
        real(r_kind) :: sigmap
        real(r_kind) :: sigman
        
        real(r_kind) :: tau
        
        ! Rematch cells on each stencil
        do iStencil = 1,nWenoStencil
          do iCell = 1,nWenoCells
            qC(iCell) = q( wenoIdx(iStencil,iCell) )
          enddo
          
          if( availableStencil(wenoType,iStencil) )then
            a(iStencil,:) = matmul( iAWENO(iStencil,:,:), qC(:) )
            if(stencil_width==3)then
              SI(iStencil) = WENO_smooth_indicator_2(a(iStencil,:))
            elseif(stencil_width>=5)then
              SI(iStencil) = WENO_smooth_indicator_3(a(iStencil,:))
            endif
          else
            a (iStencil,:) = 0
            SI(iStencil  ) = Inf
          endif
        enddo
        
        ! For WENO-Z
        tau    = 0
        iCount = 0
        do jStencil = 1,nWENOStencil-1
          if( availableStencil(wenoType,jStencil) )then
            do iStencil = jStencil+1,nWENOStencil
              if( availableStencil(wenoType,iStencil) )then
                iCount = iCount + 1
                tau = tau + abs( SI(iStencil) - SI(jStencil) )
              endif
            enddo
          endif
        enddo
        if(iCount/=0)tau = tau / iCount
        ! For WENO-Z
        
        do iPoint = 1,nWenoPoints
          if( .not.any(wenoCoef(wenoType,iPoint,:)<0) )then
            !w(:,iPoint) = wenoCoef(wenoType,iPoint,:) / ( SI + eps )**2 ! Origin
            w(:,iPoint) = wenoCoef(wenoType,iPoint,:) * ( 1. + tau / ( SI + eps ) ) ! WENO-Z
            w(:,iPoint) = w(:,iPoint) / sum(w(:,iPoint))
            
            do iTerm = 1,nWenoTerms
              wa(iPoint,iTerm) = dot_product( w(:,iPoint), a(:,iTerm) )
            enddo
            
            !qrec(iPoint) = dot_product( wa(iPoint,:), ps(iPoint,:) )
            if(present(qrec)) qrec(iPoint) = dot_product( wa(iPoint,:), ps(iPoint,:) )
            if(present(dqdx)) dqdx(iPoint) = dot_product( wa(iPoint,:), px(iPoint,:) )
            if(present(dqdy)) dqdy(iPoint) = dot_product( wa(iPoint,:), py(iPoint,:) )
          else
            ! Deal with negative linear weights, accroding to Shi J., Hu C. and Shu C.,2002
            rp = 0.5 * ( wenoCoef(wenoType,iPoint,:) + theta * abs( wenoCoef(wenoType,iPoint,:) ) )
            rn = rp - wenoCoef(wenoType,iPoint,:)
            
            sigmap = sum(rp)
            sigman = sum(rn)
            
            rp = rp / sigmap
            rn = rn / sigman
            
            !wp(:,iPoint) = rp / ( SI + eps )**2 ! Origin
            wp(:,iPoint) = rp * ( 1. + tau / ( SI + eps ) ) ! WENO-Z
            wp(:,iPoint) = wp(:,iPoint) / sum(wp(:,iPoint))
            
            !wn(:,iPoint) = rn / ( SI + eps )**2 ! Origin
            wn(:,iPoint) = rn * ( 1. + tau / ( SI + eps ) ) ! WENO-Z
            wn(:,iPoint) = wn(:,iPoint) / sum(wn(:,iPoint))
            
            do iTerm = 1,nWenoTerms
              wap(iPoint,iTerm) = dot_product( wp(:,iPoint), a(:,iTerm) )
              wan(iPoint,iTerm) = dot_product( wn(:,iPoint), a(:,iTerm) )
            enddo
            
            qrecp(iPoint) = dot_product( wap(iPoint,:), ps(iPoint,:) )
            qrecn(iPoint) = dot_product( wan(iPoint,:), ps(iPoint,:) )
            
            !qrec(iPoint) = sigmap * qrecp(iPoint) - sigman * qrecn(iPoint)
            if(present(qrec)) qrec(iPoint) = sigmap * qrecp(iPoint) - sigman * qrecn(iPoint)
            if(present(dqdx)) dqdx(iPoint) = sigmap * dot_product( wap(iPoint,:), px(iPoint,:) ) - sigman * dot_product( wan(iPoint,:), px(iPoint,:) )
            if(present(dqdy)) dqdy(iPoint) = sigmap * dot_product( wap(iPoint,:), py(iPoint,:) ) - sigman * dot_product( wan(iPoint,:), py(iPoint,:) )
          endif
        enddo
      end subroutine WENO
      
      ! 1D WENO 5th order slope limiter
      subroutine WENO5(Qrec,Q,dir,dx)
        real   (r_kind)              , intent(out) :: Qrec
        real   (r_kind), dimension(5), intent(in ) :: Q
        integer(i_kind)              , intent(in ) :: dir
        real   (r_kind), optional    , intent(in ) :: dx
        
        integer(i_kind), parameter :: nStencil = 3
        real   (r_kind), parameter :: weno_coef(3)  = [0.1, 0.6, 0.3]
        real   (r_kind), parameter :: eps           = 1.E-16
        
        real(r_kind) Qim(nStencil-1)
        real(r_kind) Qip(nStencil-1)
        real(r_kind) Qi
        
        real(r_kind), dimension(nStencil) :: stencil
        real(r_kind), dimension(nStencil) :: coefA
        real(r_kind), dimension(nStencil) :: coefB
        real(r_kind), dimension(nStencil) :: alpha
        real(r_kind), dimension(nStencil) :: beta
        real(r_kind), dimension(nStencil) :: omega
        real(r_kind), dimension(nStencil) :: xi
        real(r_kind), dimension(nStencil) :: eta
        real(r_kind), dimension(nStencil) :: etaM
        real(r_kind)                      :: xi2
        real(r_kind)                      :: lambda
        real(r_kind)                      :: lambdaM
        real(r_kind)                      :: z
        
        real(r_kind) tau40
        real(r_kind) tau41
        real(r_kind) tau5
        real(r_kind) tau
        
        integer(i_kind) iStencil
        
        Qim(2) = Q(1)
        Qim(1) = Q(2)
        Qi     = Q(3)
        Qip(1) = Q(4)
        Qip(2) = Q(5)
        
        if(dir>0)then
          stencil (1) =  Qim(2)/3. - 7./6. * Qim(1) + 11./6. * Qi     
          stencil (2) = -Qim(1)/6. + 5./6. * Qi     +  1./3. * Qip(1) 
          stencil (3) =  Qi    /3. + 5./6. * Qip(1) -  1./6. * Qip(2)
          
          coefA(1) = Qim(2) - 2. * Qim(1) + Qi
          coefA(2) = Qim(1) - 2. * Qi     + Qip(1)
          coefA(3) = Qi     - 2. * Qip(1) + Qip(2)
          
          coefB(1) =      Qim(2) - 4. * Qim(1) + 3. * Qi
          coefB(2) =      Qim(1) -      Qip(1)
          coefB(3) = 3. * Qi     - 4. * Qip(1) +      Qip(2)
        elseif(dir<0)then
          stencil (1) =  Qip(2)/3. - 7./6. * Qip(1) + 11./6. * Qi     
          stencil (2) = -Qip(1)/6. + 5./6. * Qi     +  1./3. * Qim(1) 
          stencil (3) =  Qi    /3. + 5./6. * Qim(1) -  1./6. * Qim(2)
          
          coefA(1) = Qip(2) - 2. * Qip(1) + Qi
          coefA(2) = Qip(1) - 2. * Qi     + Qim(1)
          coefA(3) = Qi     - 2. * Qim(1) + Qim(2)
          
          coefB(1) =      Qip(2) - 4. * Qip(1) + 3. * Qi
          coefB(2) =      Qip(1) -      Qim(1)
          coefB(3) = 3. * Qi     - 4. * Qim(1) +      Qim(2)
        endif
        
        beta = coefA**2 * 13. / 12. + coefB**2 * 0.25
        if(.not.any(isnan(Q)))then
          ! WENO-Z series ( WENO-Z, WENO-Z+, WENO-Z+M ) base on Luo and Wu 2021, JCP
          tau    = abs( beta(3) - beta(1) )
          xi     = ( tau + eps ) / ( beta + eps )
          if(present(dx))then
            eta    = 1._r_kind / xi
            lambda = dx**(2./3.)
            
            !! WENO-Z+
            !alpha = weno_coef * ( 1. + xi + lambda * eta )
            
            ! WENO-Z+M
            etaM     = sqrt(eta)
            xi2      = sum( weno_coef * ( 1 + xi**2 ) )
            z        = ( 1 + minval(xi) ) / xi2
            lambdaM  = lambda * z
            alpha    = weno_coef * ( 1 + xi + lambdaM * etaM )
          else
            alpha = weno_coef * ( 1 + xi )
          endif
          
          omega = alpha / sum(alpha)
          ! WENO-Z
        else
          ! origin WENO
          alpha = weno_coef / ( eps + beta )**2
          
          ! Deal with incomplete stencil
          if(dir>0)then
            if(isnan(Q(1)))then
              stencil(1) = 0
              alpha  (1) = 0
            endif
            if(isnan(Q(5)))then
              stencil(3) = 0
              alpha  (3) = 0
            endif
          endif
          
          if(dir<0)then
            if(isnan(Q(1)))then
              stencil(3) = 0
              alpha  (3) = 0
            endif
            if(isnan(Q(5)))then
              stencil(1) = 0
              alpha  (1) = 0
            endif
          endif
          
          if(isnan(Q(2)))then
            stencil(2) = 0
            alpha  (2) = 0
          endif
          if(isnan(Q(3)))then
            print*,'Center cell is NAN in WENO5'
            stop
          endif
          if(isnan(Q(4)))then
            stencil(2) = 0
            alpha  (2) = 0
          endif
          
          omega = alpha / sum(alpha)
          ! origin WENO
        endif
        
        Qrec = dot_product( stencil, omega )
        
      end subroutine WENO5
      
      ! 1D 3th order WENO limiter, accroding to 
      ! T. Lunnet et al., 2017, Monthly Weather Review
      subroutine WENO3(Qrec,Q,dir)
        real   (r_kind)              , intent(out) :: Qrec
        real   (r_kind), dimension(3), intent(in ) :: Q
        integer(i_kind)              , intent(in ) :: dir
        
        integer(i_kind), parameter :: nStencil = 2
        real   (r_kind), parameter :: weno_coef(2)  = [1./3., 2./3.]
        real   (r_kind), parameter :: eps           = 1.E-16
        
        real(r_kind) Qim(nStencil-1)
        real(r_kind) Qip(nStencil-1)
        real(r_kind) Qi
        
        real(r_kind), dimension(nStencil) :: stencil
        real(r_kind), dimension(nStencil) :: coefA
        real(r_kind), dimension(nStencil) :: coefB
        real(r_kind), dimension(nStencil) :: alpha
        real(r_kind), dimension(nStencil) :: beta
        real(r_kind), dimension(nStencil) :: omega
        
        real(r_kind) :: tau
        
        integer(i_kind) iStencil
        
        Qim(1) = Q(1)
        Qi     = Q(2)
        Qip(1) = Q(3)
        
        if(dir>0)then
          stencil (1) = -Qim(1)/2. + 3./2. * Qi
          stencil (2) =  Qi / 2. + Qip(1) / 2.
          
          coefA(1) =-Qim(1) + Qi
          coefA(2) = Qip(1) - Qi
        elseif(dir<0)then
          stencil (1) = -Qip(1)/2. + 3./2. * Qi
          stencil (2) = Qi / 2. + Qim(1) / 2.
          
          coefA(1) =-Qip(1) + Qi
          coefA(2) = Qim(1) - Qi
        endif
        
        beta = coefA**2
        
        !! origin WENO
        !alpha = weno_coef / ( eps + beta )**2
        !! origin WENO
        
        ! WENO-Z
        tau = abs( beta(1) -beta(2) )
        alpha = weno_coef * ( 1. + tau / ( eps + beta ) )
        ! WENO-Z
        
        omega = alpha / sum(alpha)
        
        Qrec = dot_product( stencil, omega )
        
      end subroutine WENO3
      
      ! Reconstruct on the left most boundary
      function left_side_recon3(q)
        real(r_kind) :: left_side_recon3
        real(r_kind),dimension(3),intent(in) :: q
        
        real(r_kind) Qip(2)
        real(r_kind) Qi
        
        Qi     = Q(1)
        Qip(1) = Q(2)
        Qip(2) = Q(3)
        
        left_side_recon3 = Qip(2)/3. - 7./6. * Qip(1) + 11./6. * Qi
        !left_side_recon3 = 1.5 * Qi - 0.5 * Qip(1)
        !left_side_recon3 = 3./8. * Qip(2) - 5./4. * Qip(1) + 15./8. * Qi
      
      end function left_side_recon3
      
      ! Reconstruct on the right most boundary
      function right_side_recon3(q)
        real(r_kind) :: right_side_recon3
        real(r_kind),dimension(3),intent(in) :: q
        
        real(r_kind) Qim(2)
        real(r_kind) Qi
        
        Qim(2) = Q(1)
        Qim(1) = Q(2)
        Qi     = Q(3)
        
        right_side_recon3 = Qim(2)/3. - 7./6. * Qim(1) + 11./6. * Qi
        !right_side_recon3 = 1.5 * Qi - 0.5 * Qim(1)
        !right_side_recon3 = 3./8. * Qim(2) - 5./4. * Qim(1) + 15./8. * Qi
      
      end function right_side_recon3
      
      function WENO_smooth_indicator_2(a) !(a_in)
        real(r_kind) :: WENO_smooth_indicator_2
        real(r_kind) :: a(:)
        !real(r_kind) :: a_in(:)
        !real(r16) :: a(lbound(a_in,1):ubound(a_in,1))
        
        !a = a_in
        
        WENO_smooth_indicator_2 = a(2)**2 + a(3)**2 + 7._r16 * a(4)**2 / 6._r16
        
      end function WENO_smooth_indicator_2
      
      function WENO_smooth_indicator_3(a) !(a_in)
        real(r_kind) :: WENO_smooth_indicator_3
        real(r_kind) :: a(:)
        !real(r_kind) :: a_in(:)
        !real(r16) :: a(lbound(a_in,1):ubound(a_in,1))
        
        !a = a_in
        
        WENO_smooth_indicator_3 = ( 720  * a(2) * a(2) + 3120 * a(3) * a(3) + 720  * a(4) * a(4) + 840   * a(5) * a(5)  &
                                  + 120  * a(4) * a(6) + 3389 * a(6) * a(6) + 3120 * a(7) * a(7) + 120   * a(2) * a(8)  &
                                  + 3389 * a(8) * a(8) + 520  * a(3) * a(9) + 520  * a(7) * a(9) + 13598 * a(9) * a(9) ) / 720
      end function WENO_smooth_indicator_3
      
      function WENO_smooth_indicator_5(a_in)
        real(r_kind) :: WENO_smooth_indicator_5
        real(r_kind) :: a_in(:)
        real(r16) :: a(lbound(a_in,1):ubound(a_in,1))
        
        a = a_in
        
        WENO_smooth_indicator_5 = a(2)**2. + (13.*a(3)**2)/3. + (3129.*a(4)**2)/80. + (87617.*a(5)**2)/140. + a(6)**2. + (7.*a(7)**2)/6. + (1./6)*a(6)*a(8) + (3389.*a(8)**2)/720. + (17./30)*a(7)*a(9) + (47459.*a(9)**2)/1120. + (1./40)*a(6)*a(10) +    &
                                  (5101.*a(8)*a(10))/1120. + (54673043.*a(10)**2)/80640. + (13.*a(11)**2)/3. + (1./24)*a(4)*a(12) + (3389.*a(12)**2)/720. + (7./20)*a(5)*a(13) + (13./18)*a(11)*a(13) + (6799.*a(13)**2)/360. + (1043./160)*a(4)*a(14) +   &
                                  (73./32)*a(12)*a(14) + (22846129.*a(14)**2)/134400. + (87617./840)*a(5)*a(15) + (13./120)*a(11)*a(15) + (306967.*a(13)*a(15))/16800. + (469977913.*a(15)**2)/172800. + (1./2)*a(6)*a(16) + (1./24)*a(8)*a(16) +          &
                                  (1./160)*a(10)*a(16) + (3129.*a(16)**2)/80. + (17./30)*a(7)*a(17) + (11./80)*a(9)*a(17) + (47459.*a(17)**2)/1120. + (1./24)*a(6)*a(18) + (73./32)*a(8)*a(18) + (24721.*a(10)*a(18))/22400. +                             &
                                  (1043./160)*a(16)*a(18) + (22846129.*a(18)**2)/134400. + (11./80)*a(7)*a(19) + (114997.*a(9)*a(19))/5600. + (114997.*a(17)*a(19))/5600. + (19583517.*a(19)**2)/12800. + (1./160)*a(6)*a(20) +                            &
                                  (24721.*a(8)*a(20))/22400. + (37850569.*a(10)*a(20))/115200. + (3129.*a(16)*a(20))/3200. + (2105043.*a(18)*a(20))/12800. + (368483712607.*a(20)**2)/15052800. + (21./5)*a(11)*a(21) + (7./20)*a(13)*a(21) +              &
                                  (21./400)*a(15)*a(21) + (87617.*a(21)**2)/140. + (1./160)*a(4)*a(22) + (5101.*a(12)*a(22))/1120. + (24721.*a(14)*a(22))/22400. + (54673043.*a(22)**2)/80640. + (21./400)*a(5)*a(23) + (7./20)*a(11)*a(23) +              &
                                  (306967.*a(13)*a(23))/16800. + (7071./800)*a(15)*a(23) + (87617./840)*a(21)*a(23) + (469977913.*a(23)**2)/172800. + (3129.*a(4)*a(24))/3200. + (24721.*a(12)*a(24))/22400. + (2105043.*a(14)*a(24))/12800. +             &
                                  (37850569.*a(22)*a(24))/115200. + (368483712607.*a(24)**2)/15052800. + (1./480)*a(2)*(240.*a(4) + 80.*a(12) + 20.*a(14) + 12.*a(22) + 3.*a(24)) + (87617.*a(5)*a(25))/5600. + (21./400)*a(11)*a(25) +                    &
                                  (7071./800)*a(13)*a(25) + (2475532433.*a(15)*a(25))/940800. + (87617.*a(21)*a(25))/5600. + (2475532433.*a(23)*a(25))/940800. + (2210903809027.*a(25)**2)/5644800. +                                                      &
                                  (a(3)*(15120.*a(5) + 2600.*a(13) + 1260.*a(15) + 390.*a(23) + 189.*a(25)))/3600
                              
      end function WENO_smooth_indicator_5
      
      function WENO_smooth_indicator_7(a_in)
        real(r_kind) :: WENO_smooth_indicator_7
        real(r_kind) :: a_in(:)
        real(r16) :: a(lbound(a_in,1):ubound(a_in,1))
        
        a = a_in
        
        WENO_smooth_indicator_7 = a(2)**2. + (13.*a(3)**2)/3. + (3129.*a(4)**2)/80. + (87617.*a(5)**2)/140. + (14127./224)*a(4)*a(6) + (252337135.*a(6)**2)/16128. + (508579./336)*a(5)*a(7) + (11102834003.*a(7)**2)/19712. +  &
                                  a(8)**2. + (7.*a(9)**2)/6. + (1./6)*a(8)*a(10) + (3389.*a(10)**2)/720. + (17./30)*a(9)*a(11) + (47459.*a(11)**2)/1120. + (1./40)*a(8)*a(12) + (5101.*a(10)*a(12))/1120. +                     &
                                  (54673043.*a(12)**2)/80640. + (47./336)*a(9)*a(13) + (34435./504)*a(11)*a(13) + (18042105247.*a(13)**2)/1064448. + (1./224)*a(8)*a(14) + (13579.*a(10)*a(14))/8064. +                         &
                                  (581814439.*a(12)*a(14))/354816. + (7505515786259.*a(14)**2)/12300288. + (13.*a(15)**2)/3. + (1./24)*a(4)*a(16) + (1./96)*a(6)*a(16) + (3389.*a(16)**2)/720. + (7./20)*a(5)*a(17) +           &
                                  (29./224)*a(7)*a(17) + (13./18)*a(15)*a(17) + (6799.*a(17)**2)/360. + (1043./160)*a(4)*a(18) + (4709./896)*a(6)*a(18) + (73./32)*a(16)*a(18) + (22846129.*a(18)**2)/134400. +                 &
                                  (87617./840)*a(5)*a(19) + (508579.*a(7)*a(19))/4032. + (13./120)*a(15)*a(19) + (306967.*a(17)*a(19))/16800. + (469977913.*a(19)**2)/172800. + (4709./896)*a(4)*a(20) +                        &
                                  (252337135.*a(6)*a(20))/96768. + (7561.*a(16)*a(20))/13440. + (18944567.*a(18)*a(20))/69120. + (82716113321.*a(20)**2)/1216512. + (508579.*a(5)*a(21))/4032. +                                &
                                  (11102834003.*a(7)*a(21))/118272. + (13./672)*a(15)*a(21) + (29183.*a(17)*a(21))/4320. + (1667122157.*a(19)*a(21))/253440. + (827161133251.*a(21)**2)/337920. +                               &
                                  (1./2)*a(8)*a(22) + (1./24)*a(10)*a(22) + (1./160)*a(12)*a(22) + (1./896)*a(14)*a(22) + (3129.*a(22)**2)/80. + (17./30)*a(9)*a(23) + (11./80)*a(11)*a(23) +                                   &
                                  (19./560)*a(13)*a(23) + (47459.*a(23)**2)/1120. + (1./24)*a(8)*a(24) + (73./32)*a(10)*a(24) + (24721.*a(12)*a(24))/22400. + (9401.*a(14)*a(24))/23040. + (1043./160)*a(22)*a(24) +            &
                                  (22846129.*a(24)**2)/134400. + (11./80)*a(9)*a(25) + (114997.*a(11)*a(25))/5600. + (190717.*a(13)*a(25))/11520. + (114997.*a(23)*a(25))/5600. + (19583517.*a(25)**2)/12800. +                 &
                                  (1./160)*a(8)*a(26) + (24721.*a(10)*a(26))/22400. + (37850569.*a(12)*a(26))/115200. + (44754957.*a(14)*a(26))/112640. + (3129.*a(22)*a(26))/3200. +                                           &
                                  (2105043.*a(24)*a(26))/12800. + (368483712607.*a(26)**2)/15052800. + (19./560)*a(9)*a(27) + (190717.*a(11)*a(27))/11520. + (138785425.*a(13)*a(27))/16896. +                                  &
                                  (45373.*a(23)*a(27))/8960. + (232084199.*a(25)*a(27))/94080. + (121599625855067.*a(27)**2)/198696960. + (1./896)*a(8)*a(28) + (9401.*a(10)*a(28))/23040. +                                    &
                                  (44754957.*a(12)*a(28))/112640. + (1732042104523.*a(14)*a(28))/5857280. + (447.*a(22)*a(28))/2560. + (18304651.*a(24)*a(28))/301056. + (3921295089347.*a(26)*a(28))/66232320. +               &
                                  (50585444357410783.*a(28)**2)/2296053760. + (21./5)*a(15)*a(29) + (7./20)*a(17)*a(29) + (21./400)*a(19)*a(29) + (3./320)*a(21)*a(29) + (87617.*a(29)**2)/140. +                               &
                                  (1./160)*a(4)*a(30) + (1./640)*a(6)*a(30) + (5101.*a(16)*a(30))/1120. + (24721.*a(18)*a(30))/22400. + (4877.*a(20)*a(30))/17920. + (54673043.*a(30)**2)/80640. +                              &
                                  (21./400)*a(5)*a(31) + (87.*a(7)*a(31))/4480. + (7./20)*a(15)*a(31) + (306967.*a(17)*a(31))/16800. + (7071./800)*a(19)*a(31) + (245947.*a(21)*a(31))/75264. +                                 &
                                  (87617./840)*a(29)*a(31) + (469977913.*a(31)**2)/172800. + (3129.*a(4)*a(32))/3200. + (14127.*a(6)*a(32))/17920. + (24721.*a(16)*a(32))/22400. + (2105043.*a(18)*a(32))/12800. +              &
                                  (199574873.*a(20)*a(32))/1505280. + (37850569.*a(30)*a(32))/115200. + (368483712607.*a(32)**2)/15052800. + (87617.*a(5)*a(33))/5600. + (508579.*a(7)*a(33))/26880. +                          &
                                  (21./400)*a(15)*a(33) + (7071./800)*a(17)*a(33) + (2475532433.*a(19)*a(33))/940800. + (6585971677.*a(21)*a(33))/2069760. + (87617.*a(29)*a(33))/5600. +                                       &
                                  (2475532433.*a(31)*a(33))/940800. + (2210903809027.*a(33)**2)/5644800. + (14127.*a(4)*a(34))/17920. + (50467427.*a(6)*a(34))/129024. + (4877.*a(16)*a(34))/17920. +                           &
                                  (199574873.*a(18)*a(34))/1505280. + (13070811329953.*a(20)*a(34))/198696960. + (365888837.*a(30)*a(34))/4515840. + (712963154333.*a(32)*a(34))/18063360. +                                    &
                                  (140082866134879519.*a(34)**2)/14306181120. + (508579.*a(5)*a(35))/26880. + (11102834003.*a(7)*a(35))/788480. + (3./320)*a(15)*a(35) + (245947.*a(17)*a(35))/75264. +                         &
                                  (6585971677.*a(19)*a(35))/2069760. + (679682189184289.*a(21)*a(35))/287006720. + (87617.*a(29)*a(35))/31360. + (878623885.*a(31)*a(35))/903168. +                                             &
                                  (282333442237789.*a(33)*a(35))/298045440. + (7284309039256637977.*a(35)**2)/20664483840. + (1./8)*a(8)*a(36) + (1./96)*a(10)*a(36) + (1./640)*a(12)*a(36) +                                   &
                                  (a(14)*a(36))/3584. + (14127./224)*a(22)*a(36) + (4709./896)*a(24)*a(36) + (14127.*a(26)*a(36))/17920. + (14127.*a(28)*a(36))/100352. + (252337135.*a(36)**2)/16128. +                        &
                                  (47./336)*a(9)*a(37) + (19./560)*a(11)*a(37) + (15.*a(13)*a(37))/1792. + (34435./504)*a(23)*a(37) + (190717.*a(25)*a(37))/11520. + (921799.*a(27)*a(37))/225792. +                            &
                                  (18042105247.*a(37)**2)/1064448. + (1./96)*a(8)*a(38) + (7561.*a(10)*a(38))/13440. + (4877.*a(12)*a(38))/17920. + (90877.*a(14)*a(38))/903168. + (4709./896)*a(22)*a(38) +                    &
                                  (18944567.*a(24)*a(38))/69120. + (199574873.*a(26)*a(38))/1505280. + (59028127.*a(28)*a(38))/1204224. + (252337135.*a(36)*a(38))/96768. + (82716113321.*a(38)**2)/1216512. +                  &
                                  (19./560)*a(9)*a(39) + (45373.*a(11)*a(39))/8960. + (921799.*a(13)*a(39))/225792. + (190717.*a(23)*a(39))/11520. + (232084199.*a(25)*a(39))/94080. +                                          &
                                  (1197465737.*a(27)*a(39))/602112. + (138785425.*a(37)*a(39))/16896. + (121599625855067.*a(39)**2)/198696960. + (1./640)*a(8)*a(40) + (4877.*a(10)*a(40))/17920. +                             &
                                  (365888837.*a(12)*a(40))/4515840. + (1297893755.*a(14)*a(40))/13246464. + (14127.*a(22)*a(40))/17920. + (199574873.*a(24)*a(40))/1505280. +                                                   &
                                  (712963154333.*a(26)*a(40))/18063360. + (22761431403211.*a(28)*a(40))/476872704. + (50467427.*a(36)*a(40))/129024. + (13070811329953.*a(38)*a(40))/198696960. +                               &
                                  (140082866134879519.*a(40)**2)/14306181120. + (15.*a(9)*a(41))/1792. + (921799.*a(11)*a(41))/225792. + (40247773253.*a(13)*a(41))/19869696. + (921799.*a(23)*a(41))/225792. +                 &
                                  (1197465737.*a(25)*a(41))/602112. + (176458381700353.*a(27)*a(41))/178827264. + (40247773253.*a(37)*a(41))/19869696. + (176458381700353.*a(39)*a(41))/178827264. +                            &
                                  (1400828669297420455.*a(41)**2)/5722472448. + (a(8)*a(42))/3584. + (90877.*a(10)*a(42))/903168. + (1297893755.*a(12)*a(42))/13246464. +                                                       &
                                  (16743073677063.*a(14)*a(42))/229605376. + (14127.*a(22)*a(42))/100352. + (59028127.*a(24)*a(42))/1204224. + (22761431403211.*a(26)*a(42))/476872704. +                                       &
                                  (293626747159200335.*a(28)*a(42))/8265793536. + (252337135.*a(36)*a(42))/3612672. + (34793506056791.*a(38)*a(42))/1430618112. + (45173351014373731.*a(40)*a(42))/1907490816. +                &
                                  (6410191990918725761813.*a(42)**2)/727389831168. + (87./56)*a(15)*a(43) + (29./224)*a(17)*a(43) + (87.*a(19)*a(43))/4480. + (87.*a(21)*a(43))/25088. +                                        &
                                  (508579./336)*a(29)*a(43) + (508579.*a(31)*a(43))/4032. + (508579.*a(33)*a(43))/26880. + (508579.*a(35)*a(43))/150528. + (11102834003.*a(43)**2)/19712. + (1./896)*a(4)*a(44) +               &
                                  (a(6)*a(44))/3584. + (13579.*a(16)*a(44))/8064. + (9401.*a(18)*a(44))/23040. + (90877.*a(20)*a(44))/903168. + (581814439.*a(30)*a(44))/354816. + (44754957.*a(32)*a(44))/112640. +            &
                                  (1297893755.*a(34)*a(44))/13246464. + (7505515786259.*a(44)**2)/12300288. + (3./320)*a(5)*a(45) + (87.*a(7)*a(45))/25088. + (29./224)*a(15)*a(45) + (29183.*a(17)*a(45))/4320. +              &
                                  (245947.*a(19)*a(45))/75264. + (181859.*a(21)*a(45))/150528. + (508579.*a(29)*a(45))/4032. + (1667122157.*a(31)*a(45))/253440. + (6585971677.*a(33)*a(45))/2069760. +                         &
                                  (140250847273.*a(35)*a(45))/119218176. + (11102834003.*a(43)*a(45))/118272. + (827161133251.*a(45)**2)/337920. + (447.*a(4)*a(46))/2560. + (14127.*a(6)*a(46))/100352. +                      &
                                  (9401.*a(16)*a(46))/23040. + (18304651.*a(18)*a(46))/301056. + (59028127.*a(20)*a(46))/1204224. + (44754957.*a(30)*a(46))/112640. + (3921295089347.*a(32)*a(46))/66232320. +                  &
                                  (22761431403211.*a(34)*a(46))/476872704. + (1732042104523.*a(44)*a(46))/5857280. + (50585444357410783.*a(46)**2)/2296053760. + (87617.*a(5)*a(47))/31360. +                                   &
                                  (508579.*a(7)*a(47))/150528. + (87.*a(15)*a(47))/4480. + (245947.*a(17)*a(47))/75264. + (878623885.*a(19)*a(47))/903168. + (140250847273.*a(21)*a(47))/119218176. +                           &
                                  (508579.*a(29)*a(47))/26880. + (6585971677.*a(31)*a(47))/2069760. + (282333442237789.*a(33)*a(47))/298045440. + (45522868146575.*a(35)*a(47))/39739392. +                                     &
                                  (11102834003.*a(43)*a(47))/788480. + (679682189184289.*a(45)*a(47))/287006720. + (7284309039256637977.*a(47)**2)/20664483840. + (14127.*a(4)*a(48))/100352. +                                 &
                                  (252337135.*a(6)*a(48))/3612672. + (90877.*a(16)*a(48))/903168. + (59028127.*a(18)*a(48))/1204224. + (34793506056791.*a(20)*a(48))/1430618112. +                                              &
                                  (1297893755.*a(30)*a(48))/13246464. + (22761431403211.*a(32)*a(48))/476872704. + (45173351014373731.*a(34)*a(48))/1907490816. + (16743073677063.*a(44)*a(48))/229605376. +                    &
                                  (293626747159200335.*a(46)*a(48))/8265793536. + (6410191990918725761813.*a(48)**2)/727389831168. +                                                                                            &
                                  (a(2)*(26880.*a(4) + 6720.*a(6) + 8960.*a(16) + 2240.*a(18) + 560.*a(20) + 1344.*a(30) + 336.*a(32) + 84.*a(34) + 240.*a(44) + 60.*a(46) + 15.*a(48)))/53760. +                               &
                                  a(3)*((21.*a(5))/5. + (87.*a(7))/56. + (13.*a(17))/18. + (7.*a(19))/20. + (29.*a(21))/224. + (13.*a(31))/120. + (21.*a(33))/400. + (87.*a(35))/4480. + (13.*a(45))/672. + (3.*a(47))/320. +   &
                                    (87.*a(49))/25088) + (508579.*a(5)*a(49))/150528. + (11102834003.*a(7)*a(49))/4415488. + (87.*a(15)*a(49))/25088. + (181859.*a(17)*a(49))/150528. +                                         &
                                  (140250847273.*a(19)*a(49))/119218176. + (452315578754789.*a(21)*a(49))/516612096. + (508579.*a(29)*a(49))/150528. + (140250847273.*a(31)*a(49))/119218176. +                                 &
                                  (45522868146575.*a(33)*a(49))/39739392. + (25839156781083324157.*a(35)*a(49))/30307909632. + (11102834003.*a(43)*a(49))/4415488. + (452315578754789.*a(45)*a(49))/516612096. +                &
                                  (25839156781083324157.*a(47)*a(49))/30307909632. + (12820383982264910635167.*a(49)**2)/40410546176._r16
      end function WENO_smooth_indicator_7
      
      subroutine get_WENO_jab(dqLdq,dqRdq,q)
        real(r_kind), dimension(5), intent(out) :: dqLdq
        real(r_kind), dimension(5), intent(out) :: dqRdq
        real(r_kind), dimension(5), intent(in ) :: q
        
        real(r_kind) :: q1
        real(r_kind) :: q2
        real(r_kind) :: q3
        real(r_kind) :: q4
        real(r_kind) :: q5
        
        real(r_kind), parameter :: C1 = 0.1_r_kind
        real(r_kind), parameter :: C2 = 0.6_r_kind
        real(r_kind), parameter :: C3 = 0.3_r_kind
        
        real(r_kind), parameter :: epsilon = 1.e-16
        
        real(r_kind) :: t7 , t20, t41, t60, t79 ,  t101,  t122,  t142, t162, &
                        t2 , t21, t42, t62, t80 ,  t110,  t123,  t143, t164, &
                        t13, t22, t43, t63, t81 ,  t102,  t124,  t144, t163, &
                        t3 , t25, t44, t88, t82 ,  t103,  t125,  t145, t165, &
                        t5 , t38, t45, t64, t85 ,  t104,  t128,  t146, t166, &
                        t10, t26, t46, t65, t86 ,  t105,  t126,  t157, t167, &
                        t4 , t27, t47, t66, t89 ,  t112,  t127,  t147, t168, &
                        t15, t28, t48, t68, t87 ,  t106,  t129,  t148, t169, &
                        t6 , t49, t50, t67, t90 ,  t107,  t130,  t149, t170, &
                        t8 , t29, t56, t69, t91 ,  t108,  t131,  t151,       &
                        t9 , t30, t57, t70, t92 ,  t111,  t132,  t150,       &
                        t11, t31, t51, t72, t93 ,  t113,  t134,  t152,       &
                        t12, t32, t52, t71, t94 ,  t115,  t133,  t160,       &
                        t14, t33, t53, t74, t95 ,  t114,  t135,  t153,       &
                        t16, t34, t73, t83, t96 ,  t116,  t136,  t154,       &
                        t23, t35, t54, t75, t97 ,  t118,  t137,  t156,       &
                        t24, t36, t55, t76, t109,  t117,  t138,  t155,       &
                        t17, t37, t58, t77, t98 ,  t119,  t139,  t158,       &
                        t18, t39, t59, t84, t99 ,  t120,  t140,  t161,       &
                        t19, t40, t61, t78, t100,  t121,  t141,  t159
        
        q1 = q(1)
        q2 = q(2)
        q3 = q(3)
        q4 = q(4)
        q5 = q(5)
        
        t7   = q2*2.0
        t2   = q1+q3-t7
        t13  = q4*2.0
        t3   = q3+q5-t13
        t5   = q3*3.0
        t10  = q2*4.0
        t4   = q1+t5-t10
        t15  = q4*4.0
        t6   = q5+t5-t15
        t8   = t2*t2
        t9   = t8*(1.3E1/1.2E1)
        t11  = t4*t4
        t12  = t11*(1.0/4.0)
        t14  = t3*t3
        t16  = t6*t6
        t23  = t14*(1.3E1/1.2E1)
        t24  = t16*(1.0/4.0)
        t17  = t9+t12-t23-t24
        t18  = abs(t17)
        t19  = epsilon+t9+t12
        t20  = 1.0/t19
        t21  = t18*t20
        t22  = t21+1.0
        t25  = q2-q4
        t38  = q3*2.0
        t26  = q2+q4-t38
        t27  = q1*(8.0/3.0)
        t28  = q3*(1.1E1/3.0)
        t49  = q2*(1.9E1/3.0)
        t29  = t27+t28-t49
        t30  = C3*t22
        t31  = epsilon+t23+t24
        t32  = 1.0/t31
        t33  = t18*t32
        t34  = t33+1.0
        t35  = C1*t34
        t36  = t25*t25
        t37  = t36*(1.0/4.0)
        t39  = t26*t26
        t40  = t39*(1.3E1/1.2E1)
        t41  = epsilon+t37+t40
        t42  = 1.0/t41
        t43  = t18*t42
        t44  = t43+1.0
        t45  = C2*t44
        t46  = t30+t35+t45
        t47  = 1.0/t46
        t48  = 1.0/(t19*t19)
        t50  = sign(1._r_kind,t17)
        t56  = t18*t29*t48
        t57  = t20*t29*t50
        t51  = t56-t57
        t52  = q2*(5.0/6.0)
        t53  = q3*(1.0/3.0)
        t73  = q1*(1.0/6.0)
        t54  = t52+t53-t73
        t55  = 1.0/(t46*t46)
        t58  = C1*t29*t32*t50
        t59  = C2*t29*t42*t50
        t61  = C3*t51
        t60  = t58+t59-t61
        t62  = q3*(1.1E1/6.0)
        t63  = q5*(1.0/3.0)
        t88  = q4*(7.0/6.0)
        t64  = t62+t63-t88
        t65  = q2*(1.0/3.0)
        t66  = q3*(5.0/6.0)
        t68  = q4*(1.0/6.0)
        t67  = t65+t66-t68
        t69  = q1*(1.9E1/3.0)
        t70  = q3*(3.1E1/3.0)
        t72  = q2*(5.0E1/3.0)
        t71  = t69+t70-t72
        t74  = t18*t48*t71
        t83  = t20*t50*t71
        t75  = t74-t83
        t76  = q2*(8.0/3.0)
        t77  = q4*(5.0/3.0)
        t84  = q3*(1.3E1/3.0)
        t78  = t76+t77-t84
        t79  = 1.0/(t41*t41)
        t80  = t18*t78*t79
        t81  = t42*t50*t71
        t82  = t80+t81
        t85  = C2*t82
        t86  = C1*t32*t50*t71
        t89  = C3*t75
        t87  = t85+t86-t89
        t90  = q1*(1.1E1/3.0)
        t91  = q2*(3.1E1/3.0)
        t92  = q4*(3.1E1/3.0)
        t93  = q5*(1.1E1/3.0)
        t94  = t90-t91+t92-t93
        t95  = q3*(2.0E1/3.0)
        t96  = t90-t91+t95
        t97  = t18*t48*t96
        t109 = t20*t50*t94
        t98  = t97-t109
        t99  = 1.0/(t31*t31)
        t100 = -t92+t93+t95
        t101 = t32*t50*t94
        t110 = t18*t99*t100
        t102 = t101-t110
        t103 = t42*t50*t94
        t104 = q2*(1.3E1/3.0)
        t105 = q4*(1.3E1/3.0)
        t112 = q3*(2.6E1/3.0)
        t106 = t104+t105-t112
        t107 = t18*t79*t106
        t108 = t103+t107
        t111 = C1*t102
        t113 = C2*t108
        t115 = C3*t98
        t114 = t111+t113-t115
        t116 = q5*(1.9E1/3.0)
        t118 = q4*(5.0E1/3.0)
        t117 = t70+t116-t118
        t119 = t18*t99*t117
        t120 = t32*t50*t117
        t121 = t119+t120
        t122 = q2*(5.0/3.0)
        t123 = q4*(8.0/3.0)
        t124 = -t84+t122+t123
        t125 = t42*t50*t117
        t128 = t18*t79*t124
        t126 = t125-t128
        t127 = C1*t121
        t129 = C3*t20*t50*t117
        t130 = C2*t126
        t131 = t127+t129+t130
        t132 = q5*(8.0/3.0)
        t134 = q4*(1.9E1/3.0)
        t133 = t28+t132-t134
        t135 = t18*t99*t133
        t136 = t32*t50*t133
        t137 = t135+t136
        t138 = C1*t137
        t139 = C3*t20*t50*t133
        t140 = C2*t42*t50*t133
        t141 = t138+t139+t140
        t142 = C1*t22
        t143 = C3*t34
        t144 = t45+t142+t143
        t145 = 1.0/t144
        t146 = q1*(1.0/3.0)
        t157 = q2*(7.0/6.0)
        t147 = t62+t146-t157
        t148 = 1.0/(t144*t144)
        t149 = C3*t29*t32*t50
        t151 = C1*t51
        t150 = t59+t149-t151
        t152 = q4*(5.0/6.0)
        t160 = q5*(1.0/6.0)
        t153 = t53+t152-t160
        t154 = q4*(1.0/3.0)
        t156 = q2*(1.0/6.0)
        t155 = t66+t154-t156
        t158 = C3*t32*t50*t71
        t161 = C1*t75
        t159 = t85+t158-t161
        t162 = C3*t102
        t164 = C1*t98
        t163 = t113+t162-t164
        t165 = C3*t121
        t166 = C1*t20*t50*t117
        t167 = t130+t165+t166
        t168 = C3*t137
        t169 = C1*t20*t50*t133
        t170 = t140+t168+t169
        
        dqLdq(1) = C3*t22*t47*(-1.0/6.0)-C3*t47*t51*t54-C3*t22*t54*t55*t60-C1*t34*t55*t60*t64-C2*t44*t55*t60*t67+C1*t29*t32*t47*t50*t64+C2*t29*t42*t47*t50*t67;
        dqLdq(2) = C3*t22*t47*(5.0/6.0)+C2*t44*t47*(1.0/3.0)+C3*t47*t54*t75-C2*t47*t67*t82+C3*t22*t54*t55*t87+C1*t34*t55*t64*t87+C2*t44*t55*t67*t87-C1*t32*t47*t50*t64*t71;
        dqLdq(3) = C3*t22*t47*(1.0/3.0)+C1*t34*t47*(1.1E1/6.0)+C2*t44*t47*(5.0/6.0)-C3*t47*t54*t98+C1*t47*t64*t102+C2*t47*t67*t108-C3*t22*t54*t55*t114-C1*t34*t55*t64*t114-C2*t44*t55*t67*t114;
        dqLdq(4) = C1*t34*t47*(-7.0/6.0)-C2*t44*t47*(1.0/6.0)+C1*t47*t64*t121+C2*t47*t67*t126-C3*t22*t54*t55*t131-C1*t34*t55*t64*t131-C2*t44*t55*t67*t131+C3*t20*t47*t50*t54*t117;
        dqLdq(5) = C1*t34*t47*(1.0/3.0)-C1*t47*t64*t137+C3*t22*t54*t55*t141+C1*t34*t55*t64*t141+C2*t44*t55*t67*t141-C3*t20*t47*t50*t54*t133-C2*t42*t47*t50*t67*t133;
        
        dqRdq(1) = C1*t22*t145*(1.0/3.0)-C1*t51*t145*t147-C1*t22*t147*t148*t150-C3*t34*t148*t150*t153-C2*t44*t148*t150*t155+C3*t29*t32*t50*t145*t153+C2*t29*t42*t50*t145*t155;
        dqRdq(2) = C1*t22*t145*(-7.0/6.0)-C2*t44*t145*(1.0/6.0)+C1*t75*t145*t147-C2*t82*t145*t155+C1*t22*t147*t148*t159+C3*t34*t148*t153*t159+C2*t44*t148*t155*t159-C3*t32*t50*t71*t145*t153;
        dqRdq(3) = C1*t22*t145*(1.1E1/6.0)+C3*t34*t145*(1.0/3.0)+C2*t44*t145*(5.0/6.0)+C3*t145*t153*(t101-t110)-C1*t98*t145*t147+C2*t108*t145*t155-C1*t22*t147*t148*t163-C3*t34*t148*t153*t163-C2*t44*t148*t155*t163;
        dqRdq(4) = C3*t34*t145*(5.0/6.0)+C2*t44*t145*(1.0/3.0)+C2*t145*t155*(t125-t128)+C3*t121*t145*t153-C1*t22*t147*t148*t167-C3*t34*t148*t153*t167-C2*t44*t148*t155*t167+C1*t20*t50*t117*t145*t147;
        dqRdq(5) = C3*t34*t145*(-1.0/6.0)-C3*t137*t145*t153+C1*t22*t147*t148*t170+C3*t34*t148*t153*t170+C2*t44*t148*t155*t170-C1*t20*t50*t133*t145*t147-C2*t42*t50*t133*t145*t155;
        
      end subroutine get_WENO_jab
      
    end module recon_mod
    
